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Abstract 

Background: An organism can respond to changing environmental conditions by adjusting gene regulation and 
by forming alternative phenotypes. In nematodes, these mechanisms are coupled because many species will form 
dauer larvae, a stress-resistant and non-aging developmental stage, when exposed to unfavorable environmental 
conditions, and execute gene expression programs that have been selected for the survival of the animal in the 
wild. These dauer larvae represent an environmentally induced, homologous developmental stage across many 
nematode species, sharing conserved morphological and physiological properties. Hence it can be expected that 
some core components of the associated transcriptional program would be conserved across species, while others 
might diverge over the course of evolution. However, transcriptional and metabolic analysis of dauer development 
has been largely restricted to Caenorhabditis elegans. Here, we use a transcriptomic approach to compare the dauer 
stage in the evolutionary model system Pristionchus pacificus with the dauer stage in C. elegans. 

Results: We have employed Agilent microarrays, which represent 20,446 P. pacificus and 20,1 43 C. elegans genes to 
show an unexpected divergence in the expression profiles of these two nematodes in dauer and dauer exit 
samples. P. pacificus and C. elegans differ in the dynamics and function of genes that are differentially expressed. We 
find that only a small number of orthologous gene pairs show similar expression pattern in the dauers of the two 
species, while the non-orthologous fraction of genes is a major contributor to the active transcriptome in dauers. 
Interestingly, many of the genes acquired by horizontal gene transfer and orphan genes in P. pacificus, are 
differentially expressed suggesting that these genes are of evolutionary and functional importance. 

Conclusion: Our data set provides a catalog for future functional investigations and indicates novel insight into 
evolutionary mechanisms. We discuss the limited conservation of core developmental and transcriptional programs 
as a common aspect of animal evolution. 

Keywords: Dauer larvae, Developmental systems drift, Transcriptomics, Evolution of gene regulation, Horizontal 
gene transfer 



Background 

The primary mechanisms of an organism to respond to 
changing environmental conditions are alterations of 
gene expression profiles and the formation of dormant 
or dormant-like developmental stages. Nematodes are 
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found in great numbers and species richness in most 
ecosystems on earth. This omnipresence is attributed to 
the evolution of a potent life history strategy to respond 
to changing environments [1]. Under favorable condi- 
tions nematodes such as Caenorhabditis elegans and 
Pristionchus pacificus undergo direct development, 
which can be completed within 3-4 days. In contrast, 
under harsh environmental conditions such as food scar- 
city, high temperature and high population density, these 
nematodes form an arrested developmental stage, the 
so-called dauer larvae [2]. 
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Dauer larvae are stress-resistant and "non-aging" and 
they facilitate the survival and dispersion of the organ- 
ism [3]. Entry into the dauer stage represents the major 
life history response of nematodes to escape unfavorable 
environmental conditions [1]. This developmental switch 
is accompanied with remarkable transcriptional and 
metabolic changes as shown in several studies in C. ele- 
gans [4-8]. Since the dauer larva is an ecologically 
induced, conserved developmental stage found in many 
species of free-living nematodes, it can be expected that 
some genetic components involved in development 
might be conserved, while those facilitating ecological 
adaptation might diverge across species. 

However, the similarity and divergence patterns of 
nematode transcriptomes have never been systematically 
investigated in the context of dauer larvae. While several 
transcriptional studies in C. elegans provide a platform 
for comparative studies [5,6,8-14], little is known in 
other free-living nematodes. Here we present a compara- 
tive transcriptomic approach to characterize the expres- 
sion profiles of the dauer and dauer exit stages in C. 
elegans and P. pacificus. We measured transcriptome- 
wide expression changes with the help of the Agilent 
microarray platform, and use this data to identify the ex- 
tent and nature of similarity or differences in the dauer- 
associated transcriptomes in the two nematode species. 

C. elegans was the first metazoan to have its complete 
genome sequenced [15]. Ever since, C. elegans has been 
at the forefront of embracing novel "omics" technolo- 
gies, including transcriptomics [16,17], proteomics 
[18,19] and RNAseq [20]. Comparative genomics is a 
powerful tool to address elementary questions in evolu- 
tionary developmental biology, such as what compo- 
nents of the regulatory "tool-kit" are conserved [21]. By 
now, the genome of six additional nematodes and many 
more transcriptomic studies have been reported (for re- 
view see [22]). Comparative transcriptomics opens up 
new vistas on how regulatory "tool-kit" components are 
employed to generate different morphologies across spe- 
cies [21,23]. In nematodes, comparative transcriptomic 
approaches provide a powerful way to analyze whole 
body responses to changes in the environment or the ex- 
posure to pathogens. The dauer stage represents a whole 
body response of the nematode to changing environ- 
mental conditions and is ideally suited for comparative 
transcriptomics studies. 

We compare C. elegans to P. pacificus, an established 
model for comparative developmental biology, evolution- 
ary biology and ecology [24]. Forward and reverse genetic 
as well as transgenic techniques have been established in 
P. pacificus, its genome has been sequenced [25] and its 
proteome has been analyzed by tandem mass spectrom- 
etry [26]. The sequencing of the P. pacificus genome 
revealed many important and unexpected features, such 



as a substantially larger size and a higher number of 
predicted protein-coding genes when compared to C. 
elegans, the presence of horizontal gene transfer and 
the duplication of genes encoding enzymes involved in the 
detoxification of xenobiotics [25]. Many of these features 
have been discussed in the context of a specific association 
with scarab beetles, in which P. pacificus is found in the 
wild [27]. Specifically, P. pacificus exists exclusively as 
dauer larvae on the living beetle and only reproduces after 
the beetle 's death by feeding on microbes that develop on 
the carcass [27-29]. This specific, so-called necromenic 
lifestyle shows the unique importance of the dauer stage 
and asks for a detailed functional investigation of the asso- 
ciated transcriptional programs. 

Two developmental processes that have been studied in 
great detail at the genetic and molecular level in P. pacificus 
are the formation of the vulva, the egg-laying structure of 
nematode females and hermaphrodites [30-32] and the 
regulation of the dauer development [33,34]. Surprisingly, 
the regulation of vulva development in C. elegans and P. 
pacificus employs different signaling pathways to control 
the formation of this homologous morphological structure. 
In dauer formation, the two transcriptional regulators 
DAF-12, a nuclear hormone receptor, and DAF-16, a fork- 
head transcription factor, are well conserved between C. 
elegans and P. pacificus [33,34]. In P. pacificus, upstream 
factors consisting of insulin and TGF-beta signaling path- 
ways in C. elegans await future analysis. Similarly, it is un- 
clear whether downstream targets that are activated by 
DAF-12 and DAF-16 are also evolutionarily conserved. To 
address this question, we have taken a transcriptomic ap- 
proach to directly compare expression profiles of the dauer 
stage (0 h) and dauer-exit stage (12 h post recovery-in- 
duction) in C. elegans and P. pacificus. 

Here we show that P. pacificus has a more dynamic tran- 
scriptome during the dauer to dauer-exit transition as com- 
pared to that in C. elegans. The expression profiles from 
the two species look surprisingly different with limited 
overlap and weak correlation between orthologous genes 
that are differentially expressed during the dauer stage and 
the dauer-recovery process. Within the conserved genes, 
functionally different GO classes and protein domains are 
enriched in each profile, most striking differences being 
observed in regulation of metabolism related genes. These 
results highlight the importance of transcriptomic studies 
in revealing functional divergence in downstream effectors 
of homologous developmental processes, despite conserva- 
tion of upstream regulatory factors. 

Results 

Significant drop in RNA abundance in dauer larvae of 
C. Elegans and P. Pacificus 

The dormant dauer stage in C. elegans is associated 
with a global repression of Pol-II based transcription, 
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down to 11 - 17% of other stages [4]. We quantified 
the amount of total RNA per worm in the dauer and 
mix-stage samples and observed that dauers of C. ele- 
gans and P. pacificus contain approximately 20-fold less 
total RNA per worm as compared to the respective 
mix-stage sample. On average, -100,000 dauer larvae 
yielded the same amount of total RNA as -5,000 mixed 
stage worms. On top of this global transcriptional re- 
pression, we detected less mRNA in the total RNA of 
dauer larvae. The mRNA proportion is reduced to 
about half of the mix-stage levels in both species as 
measured by in vitro transcription (Additional file 1: 
Figure SI a). This global repression is expected to result 
in most of the genes being down regulated in a dauer 
versus mix-stage comparison. Our normalization strat- 
egy is based on differential weighing of spiked-in RNA 
probes (see Methods, Additional file 1: Figure Sib) and 
confirms this trend for both the species (Additional file 
1: Figure Sic). 

The P. Pacificus transcriptome is more dynamic in the 
dauer to dauer-exit transition 

The Agilent microarray technology enabled simultan- 
eous expression profiling of 20,143 genes in C. elegans 
and 20,446 genes in P. pacificus. Our experimental 
setup contrasted gene expression in the dauer stage 
(time point 0 hour) with the dauer-exit samples (time 
point 12 hour) via a common reference sample (mixed 
stage) for each species. This so called "common refer- 
ence" experimental design for microarray makes use of 
a species-specific pool of RNA as a common technical 
reference, and was chosen as it potentially facilitates 
an extension of this study to other time-points and 
conditions, if needed (see [35], Methods and Add- 
itional file 1: Figure S2 for details). We have chosen 
the 12-hour post induction time-point to enable com- 
parisons with published expression profiling experi- 
ments [5]. We have chosen the same 12-hour time- 
point for the P. pacificus dauer exit samples after veri- 
fying that the post-dauer development follows similar 
kinetics in both species. We base this conclusion on 
the following morphological and developmental obser- 
vations: First, most of the population (-90%) of the 
recovering worms from C. elegans as well as P. pacifi- 
cus resume pharyngeal pumping within 3 hours after 
inducing dauer exit. Second, at the 12-hour time-point 
(the stage when we collect the dauer-exit samples), no 
discernible morphological differences are found be- 
tween the recovered dauers from both species. Third, 
the recovering animals from both species enter a 
lethargus stage between 13 to 14 hours post-recovery 
[3]. Fourth, the worms that are allowed to develop fur- 
ther at 20° C reach the next moult (L4 for C. elegans, 
J4 for P. pacificus) between 22 to 23 hours after 



induction of dauer exit. Finally, recovered worms from 
both species started laying eggs between 42 to 45 
hours post-induction. Based on all these criteria, we 
consider the 12-hour post recovery induction time- 
point to be developmentally equivalent and comparable 
across the two species. 

The number of genes, which are differentially 
expressed in the dauer to dauer exit transition (Table 1), 
is much larger in P. pacificus (4942 genes, Additional file 
2: Table SI) than that seen in C. elegans (917 genes, 
Additional file 3: Table S2), suggesting more dynamic 
transcriptional changes accompanying dauer recovery in 
P. pacificus. Also, the number of up-regulated versus 
down-regulated genes in the dauers to dauer-exit com- 
parison is much larger in P. pacificus than that in C. ele- 
gans, where the number of up- and down- regulated 
genes is essentially the same (Table 1). Thus P. pacificus 
dauers potentially require the activity of more genes to 
survive in their ecological niche than the dauer recovery 
stage. 

Signatures of conservation and divergence in the 
transcriptomes of the two species 

Using pairwise best BLAST mapping, we identified 
6,126 1:1 orthologous gene pairs in the two species, 
which are represented on both the arrays, while the 
remaining 14,212 P. pacificus genes and 13,099 C. ele- 
gans genes represent unresolved homology relations 
and lineage-specific genes. Within the 6,126 orthologs, 
we find that only 184 gene pairs are expressed at sig- 
nificantly different levels in the dauer versus dauer-exit 
comparisons of both species (FDR corrected p-value < 
0.05) (Figure la). This small overlap is nonetheless 
statistically significant (Fisher's exact test p-value = 
0.029) indicating some "conservation" between dauer 
related genes in the two species. 

The Pearson's correlation coefficient (r) between the 
fold changes of the 1:1 orthologs is r = 0.12 when calcu- 
lated over all 6,126 orthologs, and increases to r = 0.29 
when calculated on the 184 genes that are significantly 

Table 1 Differentially expressed genes in the dauer 
versus dauer-exit comparison 

P. pacificus transcriptome is relatively more dynamic in dauer to 
dauer-exit transition (FDR corrected p-value < = 0.05) 



C. elegans 
P. pacificus 



Up 


Down 


Total 


Total genes on array 


476 


441 


917 


20143 


3545 


1394 


4939 


20446 



The Agilent arrays interrogate 20,446 genes for P. pacificus and 20,143 genes 
for C elegans. We measure expression changes in the dauer to dauer-exit 
(12-hour time-point) transition using species-specific mix-stage sample as a 
common technical reference. Genes were called significantly differentially 
expressed (up- or down- regulated) based on a FDR corrected p-value cut-off 
of 0.05. P. pacificus transcriptome appears to be more dynamic during this 
transition. 
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(a) Differentially expressed genes in C. elegans and P. pacificus 
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(b) C. elegans / P. pacificus, 1:1 orthologs (n= 6126 ; signif= 184 ) 
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Pearson correlation coefficient: all genes r = 0.12 ; significant genes r = 0.29 

Figure 1 Limited conservation in the dauer related transcriptomes of the two species, (a) Overlap between genes differentially expressed 
in P. pacificus (green) and C. elegans dauers versus dauer-exit (12 hours) samples. The rectangular boxes represent the entire transcriptome on the 
array and their area of overlap represents the 1:1 orthologs between the two species, which are represented on both microarrays. The ovals 
represent the set of differentially expressed genes in each species. 184 orthologous gene pairs are called significantly differentially expressed in 
both the species, indicating limited conservation (Fisher's exact test p-value = 0.029). A substantial number of non-orthologous, species-specific 
genes are also differentially expressed in both the species. These gene sets are further analyzed by dividing them into distinct classes which are 
indicated in blue text (see main text for details) (b) Comparison of log2 expression fold changes for the set of 6,1 26 1 :1 orthologs for the dauer 
versus dauer-exit comparison. Pearson's correlation coefficient over the entire set of 6,126 orthologs is r = 0.1 2, and increases to r = 0.29 for the 
184 orthologs, which are significantly differentially expressed in both species (FDR corrected p-value<=0 0.05; purple boxes mark the significant 
orthologs, number per quadrant I - IV: 51 +51 +65 + 17 = 184). 



differential in both species (Figure lb). To further 
characterize the cross-species similarities and differences 
in transcriptomes, we group the set of 1:1 orthologs into 
the following classes: 



(1) Class 1 = common orthologs that show a 
concordant pattern of differential regulation 
between the two species (up-regulated in both 
species, n = 51, quadrant 1, and down-regulated in 
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both species n = 65, quadrant 3, Figure lb. 
Total = 116 out of 184 common orthologs). 
These are expected to be the most conserved 
part of the transcriptomes across species 
(r = 0.88) 

(2) Class 2 = orthologous genes that show a negative 
correlation in their direction of fold- change 

(n = 51, quadrant 2, and n = 17, quadrant 4, 
Figure lb. Total = 68 out of total 184 common 
orthologs). This class of genes exhibits the most 
divergent expression pattern in the two 
transcriptomes (r = -0.80) and is comparable in size 
to class 1. 

(3) Class 3 = Orthologs differentially expressed only in 
C. elegans dauers (n = 178 in C. elegans, Figure la) 

(4) Class 4 = Orthologs differentially expressed only 
in P. pacificus expression profile (n = 2626 in 

P. pacificus, Figure la). 

Gene ontology term enrichment 

To better understand which functions are enriched in 
these gene sets, we made use of Gene Ontology (GO) 
term annotations from WormBase [36]. The set of 6,126 
1:1 orthologs was used as the background set, and GO 
annotations were assigned to P. pacificus genes by 



directly mapping them from the corresponding C. ele- 
gans 1:1 orthologs. Enrichment statistics were calculated 
with the Bioconductor package topGO [37] (see Meth- 
ods). We found the biological process of "neuropeptide- 
signaling pathway" to be the most significantly enriched 
(Table 2) in the class 1 gene set, which subsumes the 
most conserved expression pattern across the two spe- 
cies. We find six out of the 16 genes belonging to this 
GO class to be upregulated in dauers of both C. elegans 
and P. pacificus, and four of them are members of 
"FMRF-Like Peptide" family of neuropeptides, which 
also have a known role in regulating pharyngeal pump- 
ing [38]. GO enrichment analysis of class 2 genes that 
show an opposite fold change in their expression pattern 
between species yields "molting cycle, collagen and 
cuticulin-based cuticle" as the most significantly 
enriched term (Table 3). This suggests that the dauers of 
both species possibly use different genetic components 
for synthesis and/or shedding of the cuticle. This poten- 
tial difference in molting related cuticle processing in 
the two species is further supported by the observation 
that this GO term is also the most significantly enriched 
in Class 3 genes - the set of orthologous genes differen- 
tially regulated exclusively in C. elegans (Table 4). For 
class 4 genes, the most significantly enriched GO terms 



Table 2 GO terms enriched in orthologs with same direction of fold change in both species 

GO:BP terms enriched in orthologs regulated in same direction in both species 



GO.ID 


Term 


Annotated 


Significant Expected 


p-Value 


GO:0007218 


neuropeptide signaling pathway 


16 


6 


0.29 


2.00E-07 


GO:0055114 


oxidation-reduction process 


219 


10 


3.97 


0.0056 


GO:0006576 


cellular biogenic amine metabolism 


/ 


2 


0.13 


0.0064 


GO:0006508 


proteolysis 


197 


9 


3.57 


0.0086 


GO:0046496 


nicotinamide nucleotide metabolism 


10 


2 


0.18 


0.0133 


GO:0009069 


serine family amino acid metabolism 


10 


2 


0.18 


0.0133 


GO:0019752 


carboxylic acid metabolism 


128 


8 


2.32 


0.0172 


GO:0009309 


amine biosynthetic process 


33 


3 


0.6 


0.0212 


GO:0040019 


positive regulation of embryonic development 


13 


2 


0.24 


0.0222 


GO:0008202 


steroid metabolic process 


14 


2 


0.25 


0.0256 


GO:0008033 


tRNA processing 


19 


2 


0.34 


0.0455 


GO:0006006 


glucose metabolic process 


20 


2 


0.36 


0.0499 


GO:MF terms enriched in orthologs regulated in same direction in both species 










GO.ID 


Term 


Annotated 


Significant Expected 


p-Value 


GO:0016616 


oxidoreductase activity; acting on the CH-OH group of donors, NAD or NADP as acceptor 


34 


4 


0.64 


0.0034 


GO:00 16491 


oxidoreductase activity 


283 


1-1 


5.29 


0.0107 


GO:0016747 transferase activity; transferring acyl groups other than amino-acyl groups 


53 


4 


0.99 


0.0164 


GO:00 16831 


carboxy-lyase activity 


11 


2 


0.21 


0.017 


GO:00041 75 


endopeptidase activity 


115 


6 


2.15 


0.0193 


GO:0008233 


peptidase activity 


180 


10 


3.37 


0.0258 


GO:0051287 


NAD binding 


14 


2 


0.26 


0.0271 
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Table 3 GO terms enriched in orthologs with opposite direction of fold change in both species 


G0:BP terms enriched in orthologs regulated in opposite direction in both species 


GO.ID Term 


Annotated Significant Expected 


p-Value 


GO:0018996 molting cycle; collagen and cuticulin-based cuticle 




157 7 


1.77 


0.0016 


GO:0010171 body morphogenesis 




322 10 


3.63 


0.0025 


GO:004001 1 locomotion 




789 17 


8.89 


0.0038 


GO:0040010 positive regulation of growth rate 




1 024 18 


11.54 


0.0244 


GO:0040018 positive regulation of multicellular organism growth 




167 5 


1.88 


0.0382 


GO:0019318 hexose metabolic process 




28 2 


0.32 


0.039 


GO:MF terms enriched in orthologs regulated in opposite direction in both species 










GO.ID Term 


Annotated Significant Expected 


p-Value 


GQ0042302 structural constituent of cuticle 




21 4 


0.2 


3.80E-05 


GQ0015078 hydrogen ion transmembrane transporter activity 




30 2 


0.29 


0.033 


GQ0050662 coenzyme binding 




83 3 


0.8 


0.045 


GO:0046912 transferase activity; transferring acylgroups, acyl groups converted into alkyl on transfer 


5 1 


0.05 


0.047 


are mostly development related but are too broad to originate from other published microarray data sets. To 
point to any specific function (Table 5). this end, we retrieved predefined expression clusters from 

WormBase [36] and used them for a meta-level analysis. 
Expression cluster analysis These expression clusters represent sets of genes reported 
Apart from analyzing enrichment for GO annotations, we to be co-expressed under various conditions. We 
also performed a functional analysis based on overlap of restricted the analysis to the set of 1:1 orthologs, and 
our data with gene expression clusters in C. elegans, which assigned expression clusters to P. pacificus genes by 


Table 4 GO terms enriched in orthologs, which are only differentially expressed in C. elegans 






GO:BP terms enriched in orthologs differentially expressed exclusively in C. elegans 


GO.ID Term 


Annotated 


Significant 


Expected 


p-Value 


GO:0018996 molting cycle; collagen and cuticulin-based cuticle. 


157 


17 


4.19 


5.70E-07 


GO:004001 1 locomotion 


789 


34 


21.07 


0.0019 


GO:0040018 positive regulation of multicellular organism growth 


167 


11 


4.46 


0.0046 


GO:0006560 proline metabolic process 


5 


2 


0.13 


0.0067 


GO:0010171 body morphogenesis 


322 


16 


8.6 


0.0106 


GO:0009084 glutamine family amino acid biosynthesis 


8 


2 


0.21 


0.0178 


GO:0006694 steroid biosynthetic process 


10 


2 


0.27 


0.0277 


GO:0055114 oxidation-reduction process 


219 


11 


5.85 


0.0306 


GO:0006508 proteolysis 


197 


10 


5.26 


0.0362 


GO:0018991 oviposition 


145 


8 


3.87 


0.0387 


GO:MF terms enriched in orthologs differentially expressed exclusively in C. elegans 










GO.ID Term 


Annotated 


Significant 


Expected 


p-Value 


GQ0042302 structural constituent of cuticle 


21 


10 


0.53 


1.60E-11 


GO:0004222 metalloendopeptidase activity 


50 


6 


1.25 


0.0014 


GQ0016776 phosphotransferase activity; phosphate group as acceptor 


5 


2 


0.13 


0.0059 


GO:0016903 oxidoreductase activity; acting on the aldehyde or oxo group of donors 


8 


2 


0.2 


0.0157 


GO:0003854 3-beta-hydroxy-delta5-steroid dehydrogenase activity 


8 


2 


0.2 


0.0157 


GO:0019205 nucleobase; nucleoside; nucleotide kinase activity 


8 


2 


0.2 


0.0157 


GO:0020037 heme binding 


51 


4 


1.28 


0.0374 


GO:0004601 peroxidase activity 


14 


2 


0.35 


0.0464 
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Table 5 GO terms enriched in orthologs, which are only differentially expressed in P. pacificus 
G0:BP terms enriched in orthologs differentially expressed exclusively in P. pacificus 



GO.ID 


Term 


Annotated 


Significant 


Expected 


p-Value 


GO:0006898 


receptor-mediated endocytosis 


433 


231 


192.89 


6.30E-05 


GO:0009792 


embryo development ending in birth or egg hatching 


1535 


740 


683.81 


0.00015 


GO:0002119 


nematode larval development 


1123 


549 


500.27 


0.00034 


GO:0040020 


regulation of meiosis 


45 


30 


20.05 


0.00219 


GO:0042127 


regulation of cell proliferation 


45 


30 


20.05 


0.00219 


GO:0006412 


translation 


165 


93 


73.5 


0.00691 


GO:00 16246 


RNA interference 


63 


38 


28.07 


0.00812 


GO:0055114 


oxidation-reduction process 


219 


115 


97.56 


0.00912 


GO:0000003 


reproduction 


1235 


585 


550.17 


0.00934 


GO:0040007 


growth 


1481 


694 


659.75 


0.01358 


GO:0009396 


folic acid-containing compound biosynthetic process 


5 


5 


2.23 


0.01749 


GO:0009067 


aspartate family amino acid biosynthesis 


10 


8 


4.45 


0.02546 


GO:0009072 


aromatic amino acid family metabolic process 


10 


8 


4.45 


0.02546 


GO:0006732 


coenzyme metabolic process 


41 


27 


18.26 


0.0326 


GO:0051603 


proteolysis involved in cellularprotein catabolic process 


36 


24 


16.04 


0.03319 


GO:0006396 


RNA processing 


72 


40 


32.07 


0.03822 


GO:0006511 


ubiquitin-dependent protein catabolic process 


29 


18 


12.92 


0.04326 


GO:0006418 


tRNA aminoacylation for protein translation 


33 


20 


14.7 


0.0461 1 


GO:0015684 


ferrous iron transport 


9 


/ 


4.01 


0.04679 


GO:0042026 


protein refolding 


9 


/ 


4.01 


0.04679 


GO:0006399 


tRNA metabolic process 


51 


32 


22.72 


0.04805 


G0:MF terms enriched in orthologs differentially expressed exclusively in P. pacificus 








GO.ID 


Term 


Annotated 


Significant 


Expected 


p-Value 


GO:0004298 


threonine-type endopeptidase activity 


13 


12 


5.76 


4.30E-04 


GO:0003899 


DNA-directed RNA polymerase activity 


22 


1/ 


9.74 


0.00168 


GO:00 16491 


oxidoreductase activity 


283 


149 


125.32 


0.00192 


GO:0008026 


ATP-dependent helicase activity 


49 


32 


21.7 


0.00229 


GO:0050662 


coenzyme binding 


83 


48 


36.76 


0.00828 


GO:0003735 


structural constituent of ribosome 


93 


53 


41.18 


0.00847 


GO:0008168 


methyltransferase activity 


67 


39 


29.67 


0.01438 


GO:00 16884 


carbon-nitrogen ligase activity; with glutamine as amido-N-donor 


10 


8 


4.43 


0.02443 


GO:0051082 


unfolded protein binding 


28 


18 


12.4 


0.02591 


GO:0030170 


pyridoxal phosphate binding 


34 


21 


15.06 


0.02977 


GO:0004812 


aminoacyl-tRNA ligase activity 


34 


21 


15.06 


0.02977 


GO:0003993 


acid phosphatase activity 


12 


9 


5.31 


0.03154 


GO:0031072 


heat shock protein binding 


23 


15 


10.19 


0.03483 


GO:00041 90 


aspartic-type endopeptidase activity 


9 


/ 


3.99 


0.04517 


GO:0015093 


ferrous iron transmembrane transporter activity 


9 


/ 


3.99 


0.04517 


GO:0008483 


transaminase activity 


9 


7 


3.99 


0.04517 



mapping annotation from the corresponding C. elegans 
orthologs. The significance scores for enrichment of each 
pre-defined cluster was calculated separately for the sig- 
nificantly up- and down- regulated genes from the dauer 



versus dauer exit comparison in our data (see Methods for 
details), and are summarized in Table 6. 

These significantly enriched clusters enable us not 
only to identify system-wide trends in our data but 
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also validate them against existing data-sets. For ex- 
ample, we see a highly significant overlap between 
our C. elegans data and previously reported expression 
profiles of dauers and dauer-exit stages [5,9]. Specific- 
ally, the cluster "Dauer enriched" identified by Wang 
and Kim in their time-course study [5] is enriched in 
genes called up-regulated in our data (Table 6, row 1, 
column "cel_dauer"). Furthermore, other clusters from 
later dauer-exit time points in the same study show a 
significant overlap with genes called over-expressed in 
our dauer-exit samples (Table 6, rows 2 to 4, column 
"cel_exit"), thus indicating good agreement between 
the two data-sets. Interestingly, the clusters corre- 
sponding to "early" and "climbing" genes in C. elegans 
(Table 6, rows 2 and 3) also show a significant over- 
lap with genes over-expressed in our P. pacificus 
dauer-exit samples (column "ppa exit" in Table 6), 
pointing towards a conservation of a part of the ac- 
tive transcriptome during dauer recovery in the two 
species. A similar trend is also seen for clusters 
obtained from a microarray study of dauers from daf- 
c/TGF-beta mutants [9] (Table 6, rows 5 and 6). We 
also see a significant overlap between genes regulated 
in response to heat shock, oxidative stress, osmotic 
stress and response to pathogens, highlighting the fact 
that activation of stress response pathways is a com- 
mon feature of dauer stage in both C. elegans and P. 
pacificus. Further, genes activated in response to star- 
vation and on resumption of feeding are also found 
to be enriched in dauer and dauer-exit stages respect- 
ively for both the species, in good agreement with 
non-feeding status of dauers and resumption of a 
feeding program during dauer-exit. We observe a sig- 
nificant overlap between the expression clusters that 
are enriched in the dauers of both the species (Fish- 
er's exact test p-value = 0.0185), indicating a conserved 
signature of gene expression. This overlap of expres- 
sion clusters is stronger for the dauer-exit stages of 
both the species (Fisher's exact test p-value = 6.833E- 
06), indicating that similar transcriptional programs 
are activated during dauer recovery in both the spe- 
cies. In summary, in spite of the low overlap and cor- 
relation between differentially expressed genes in the 
two species, functional analysis based on GO and ex- 
pression cluster enrichment could identify conserved 
aspects of the transcriptional program associated with 
the dauer recovery process in the two species. 

Apart from the differentially expressed genes in the 
limited set of 1:1 orthologous genes, there is a com- 
parable number of differentially expressed genes for 
which either no orthologs exist in the other species 
or the orthology cannot be uniquely resolved (2,129 
P. pacificus specific genes versus 2626 + 184 = 2810 1:1 
orthologs in P. pacificus, 555 C. elegans specific genes 



versus 178 +184 = 362 1:1 orthologs in C. elegans 
data, Figure la). This indicates that the non- 
conserved part of the active transcriptome show a 
similar dynamic like the conserved part and might 
play a significant role in the biology of dauer larvae 
of the two species. Further functional studies would 
be required to understand the adaptive significance of 
these species-specific genes. 

PFAM protein domain based analysis identifies potentially 
conserved and diverged functional gene classes 

We further investigated the functional differences in 
the two transcriptomes by stratifying gene expression 
based on their protein domain annotations. This ap- 
proach has the benefit of not being dependent on 
orthology relationships (like some of the analyses per- 
formed above), and can identify conserved functional 
signatures if the same gene function is carried out by 
different paralogs in the two species. For this analysis, 
we inferred possible gene functions by annotating 
13,344 gene loci in C. elegans and 14,018 gene loci in 
P. pacificus with PFAM protein domains (p-value of 
domain match < 0.001). We then stratified the gene 
sets from the two species into gene families by PFAM 
domains with at least five members on each micro- 
array. These 441 gene families may partially overlap 
and are represented by their median log2 fold expres- 
sion changes. The Pearson's correlation coefficient of 
log2 fold changes is r = 0.52 for the dauer versus 
dauer exit comparison (Figure 2), an improvement 
over the aforementioned correlation values for 1:1 
orthologs (Figure lb). Some protein domains indeed 
show concordant expression in both the species, such 
as the HSP20 family of heat shock proteins that is 
induced in dauers of both the species (Figure 2) and 
potentially confers stress resistance. Similarly, the 
Peptidase_S28 family of proteins is repressed in the 
dauer stage of both species possibly because it con- 
tains proteins with lysosomal activity, which might 
not be required during the repressed metabolic state 
of dauers [7,39]. Interestingly, we also find some pro- 
tein domains with totally opposite expression pat- 
terns, such as some of the proteasome subunit 
domains. The genes with domains "Proteasome" and 
"Proteasome A N" are induced strongly in P. pacifi- 
cus dauers but are down regulated in C. elegans 
dauers. Since proteasome function has been impli- 
cated in regulation of longevity [40-42], up regulation 
of these specific proteasome domains in P. pacificus 
might be correlated with the increased longevity of P. 
pacificus dauers versus C. elegans dauers [43]. Thus, 
in summary, PFAM annotation based analysis of the two 
transcriptomes also identifies some conserved functional 
signatures as well as some divergent signatures. This 
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Table 6 Expression cluster enrichment analysis 



Expression clusters over-represented in "dauer enriched" and "exit enriched genes" in C. elegans and P. pacificus 


Expression clusters related to dauer larvae, stress response and starvation / feeding 


SI.No. 


Expression Cluster 


cel_dauer 


cel_exit 


ppa_dauer 


ppa_exit 


1 


Wang_Kim_WBPaper00005859_DauerEnriched 


10.45 


0 


0 


0 


2 


Wang_Kim_WBPaper00005859_EarlyGenes 


0 


5.99 


0 


10.67 


3 


Wang_Kim_WBPaper00005859_ClimbingGenes 


0 


14.58 


0 


4.62 


4 


Wang_Kim_WBPaper00005859_LateGenes 


0 


63.09 


0 


0 


5 


WBPaper00024393:strongly_regulated_dauer_genes_UP 


1.39 


0 


0 


0 


6 


WBPaper00024393:strongly_regulated_dauer_genes_DOWN 


0 


8.49 


0 


10.85 


/ 


WBPaper00034757:up_by_oxidative_stress 


4.2 


0 


0 


0 


8 


WBPaper00034757:down_by_oxidative_stress 


0 


6.47 


0 


6.82 


9 


WBPaper00035227:heat_shock_regulated 


9.74 


0 


0 


0 


10 


WBPaper00035873:dpy-10_regulated 


0 


25.06 


0 


3.89 


11 


WBPaper00035873:dpy-9_regulated 


0 


2.02 


0 


1.57 


12 


WBPaper00035873:osm-11_regulated 


0 


2.02 


0 


1.57 


13 


WBPaper00035873:osm-7_regulated 


0 


2.4 


0 


2.07 


14 


WBPaper00035873:osm-8_regulated 


0 


42.83 


0 


445 


15 


WBPaper00035873:osmotically_regulated 


0 


0 


0 


4.65 


16 


WBPaper00032948:StarveUp2 


6.17 


0 


2.5 


0 


1/ 


WBPaper00032948:StarveUp3 


7.15 


0 


0 


0 


18 


WBPaper00032948:StarveUp4 


10.45 


0 


0 


0 


19 


WBPaper00032948:FedUp 


0 


78.95 


0 


5.12 


20 


WBPaper00032948:MoltOssilate 


0 


74.78 


0 


0 


21 


WBPaper00032062:age_regulated_genes 


0 


5.37 


0 


4.75 


Expression clusters related to to pathogen response, RNAi machinery etc. 










SI.No. 


Expression Cluster 


cel_dauer 


cel_exit 


ppa_dauer 


ppa_exit 


22 


WBPaper00028482:PA14_upregulate 


0 


0 


0 


2.03 


23 


WBPaper00028789:PA14_vs_gacA_downregulated_4hr 


0 


5.9 


0 


0 


24 


WBPaper00028789:PA14_vs_gacA_downregulated_8hr 


0 


6.7 


0 


0 


25 


WBPaper00028789:PA14_vs_OP50_downregulated_8hr 


0 


6.85 


0 


1.81 


26 


WBPaper00028789:PA14_vs_OP50_upregulated_4hr 


0 


0 


0 


1.4 


27 


WBPaper00028789:PA14_vs_OP50_upregulated_8hr 


0 


0 


0 


2.02 


28 


WBPaper00030985:Enterococcus_faecalis_upregulated 


0 


0 


0 


1.6 


29 


WBPaper00028789:pmk-1_downregulated 


0 


0 


0 


1.46 


30 


WBPaper00029437:dcr-1_upregulated 


3.13 


0 


0 


0 


31 


WBPaper00029437:rde-4_upregulated 


0 


0 


0 


1.69 


32 


WBPaper000271 1 1 :eri-1 (mg366)_downregulated 


0 


2.23 


0 


0 


33 


WBPaper000271 1 1 :rde-3(r459)_upregulated 


0 


0 


0 


1.52 


34 


WBPaper00035892:KIM5_regulated 


0 


0 


0 


4.79 


35 


WBPaper00035892:KIM5_vs_OP50_Up 


0 


0 


0 


2.36 


Tissue specific, and other expression clusters 










SI.No. 


Expression Cluster 


cel_dauer 


cel_exit 


ppa_dauer 


ppa_exit 


36 


WBPaper00030839:Embryo_Pan_Neuronal 


0 


0 


6.67 


0 


37 


WBPaper00030839:Larval_Pan_Neuronal 


4.75 


0 


2.71 


0 


38 


WBPaper00031003:0hr_muscle_depleted 


0 


2.99 


0 


0 
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Table 6 Expression cluster enrichment analysis (Continued) 



39 


WBPaper00031003:24hr_muscle_depleted 


0 


2.11 


0 


2.73 


40 


WBPaper00031003:total_muscle_depleted 


0 


2.29 


0 


0 


41 


WBPaper00031003:total_muscle_enriched 


0 


0 


242 


0 


42 


WBPaper00031532:Larva_Pan_Neuronal_Depleted 


0 


28.03 


0 


11.41 


43 


WBPaper00031532:Larva_Pan_Neuronal_Enriched 


2.48 


0 


4.69 


0 


44 


WBPaper00026980:intestine_enriched 


0 


0 


346 


1.5 


'15 


WBPaper00024671:AFD_AWB_vs_unsorted_downregulated 


0 


3.02 


0 


3.59 


46 


WBPaper00031832:slr-2_regulated 


0 


0 


0 


2.75 


47 


WBPaper00033101:spr-5_regulated 


0 


6.44 


0 


0 


48 


WBPaper00034739:N2lessDR1 350 


0 


39.13 


0 


0 


49 


WBPaper00034739:RIL17lessRIL14 


0 


18.11 


0 


0 


50 


WBPaper00035905:FBF-1_Associated 


0 


0 


5.54 


0 


51 


WBPaper0003761 1 :RNP-8-associated 


0 


0 


4.92 


0 


52 


WBPaper00025032:PAL-1_target_genes 


0 


2.58 


0 


0 



This table lists - loglO transformed p-values of tests for enrichment of dauer and dauer-exit genes from both species for 52 out of 169 expression clusters. These 
52 expression clusters passed the FDR corrected p-value threshold of 0.05 for at least one condition. The overlap between dauer enriched clusters and dauer-exit 
enriched clusters is statistically significant (Fisher's exact test p-values = 0.01 85 and 6.833E-06 respectively). 



analysis also provides an important catalog of gene func- 
tion, and the role of genes corresponding to the protein 
families identified here can be studied in more detail in P. 
pacificus dauers in the future. 



Metabolic recovery during dauer exit displays different 
expression patterns of key enzymes in P. Pacificus 

Dauer larvae of C. elegans are known to undergo a re- 
markable shift in their intermediary metabolism to 



C. elegans / P. pacificus Microarray 1 :1 PFAM families (n= 441 ) 



Proteasome ■ 



Proteasome_A_N 



PCI 



■ Peptidase_S28 



HSP20 



4 




-2.0 



-1.5 



-1.0 



-0.5 



0.0 



0.5 



1.0 



C. elegans Dauer/DauerExit log2 FC 
Pearson correlation coefficient between significant PFAM domains r = 0.47 

Figure 2 Log2 fold changes of median PFAM gene expression levels in P. pacificus and C. elegans. All genes on the microarrays were 
stratified by PFAM domain annotation into partially overlapping gene families with at least 5 members. A positive log2 fold change indicates that 
the corresponding genes were up-regulated in the dauer stages in comparison to the dauer exit time point. 
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survive long periods without actively feeding [7,39]. We 
therefore looked in more detail into expression changes 
in metabolism related genes using the KEGG database 
[44]. Gene-pathway assignments were obtained for C. 
elegans from KEGG and transferred to P. pacificus via 
the 6,126 orthologs identified through best pairwise 
BLAST mapping. For each KEGG pathway, we deter- 
mined the number of genes from the 1:1 orthologous set 
that are upregulated, downregulated or show no differ- 
ences in expression levels in the dauer versus dauer-exit 
comparisons in both species. We observed that the three 
key processes of central carbon metabolism, namely gly- 
colysis, the Krebs (TCA) cycle and oxidative phosphoryl- 
ation, show entirely different dynamics during the 
dauer-exit time course in the two species (see Figure 3). 
For C. elegans, mRNA abundance levels for most of the 
genes in the three pathways are found to be similar in 
the dauers and dauer-exit stage (category "no_change", 
Figure 3). Only a few genes are downregulated, while 
none of the genes is upregulated in C.elegans. In con- 
trast, we see a more dynamic regulation of metabolic 
pathways in P. pacificus, with comparatively larger 



number of downregulated genes as well as a few upregu- 
lated genes. These differences in metabolic recovery in 
P. pacificus could have an obvious and trivial explan- 
ation in that P. pacificus might have a different develop- 
mental rate compared to C. elegans. However, as 
described above, several morphological and developmen- 
tal traits suggest that the course of post-dauer develop- 
ment in P. pacificus and C. elegans is very similar. 
Hence, we hypothesize that metabolic regulation in P. 
pacificus dauers might itself be inherently different from 
that in C. elegans. This would also be consistent with 
the dauer recovery studies showing that P. pacificus 
dauers can survive much longer than C. elegans dauers 
[43], one of the possible reason being their ability to 
metabolize stored fats at a different rate and/or in a dif- 
ferent way. Our genome-wide expression studies reflect 
this potential difference as divergence in pathway ex- 
pression profiles between the two species. Since gene ex- 
pression profiles are dynamic and can be sensitive to 
differences in developmental timing (e.g. [11]), more 
fine-grained and detailed functional studies will be 
needed to verify whether these differences are really due 



C. elegans 



P. pacificus 



Glucose 



Glycolysis / Gluconeogenesis 


C. elegans 


up 


down 


no_change 


total 


0 


4 


21 


25 



Glycolysis / Gluconeogenesis 


P. pacificus 


up 


down 


no_change 


total 


3 


6 


13 


22 



Pyruvate -►Acetyl-CoA 



Citrate cycle 


C. elegans 


up 


down 


no_change 


total 


0 


6 


17 


23' 4 










* 






















Oxidative phosphorylation 


C. elegans 


up 


down 


no_change 


total 


0 


3 


70 


73 



FADH 
NADH 

0 2 



I 



Citrate cycle 


P. pacificus 


up 


down 


no_change 


total 


,' 3 


12 


8 


23 


t 










t 










Oxidative phosphorylation 


P. pacificus 


up 


down 


no_change 


total 


17 


10 


34 


61 



ATP 

Figure 3 Differences in regulation of metabolism related genes during the dauer exit time course in C. elegans and P. pacificus. 

Gene-to-pathway assignments were retrieved for C. elegans from the KEGG database. Pathway assignments were transferred to P. pacificus by 
mapping 1:1 orthologs. All tables show genes that are up, down or at a similar expression level in a dauer vs. dauer exit (12 hours post induction) 
comparison. 



Sinha et al. BMC Genomics 2012, 13:254 
http://www.biomedcentral.com/1471 -21 64/1 3/254 



Page 12 of 17 



to inherent metabolic differences or an artifact of slight 
variations in developmental timing. 

Genes acquired by horizontal gene transfer and orphan 
genes in P. Pacificus are developmental^ regulated 

Previous studies have identified genes in P. pacificus that 
likely originate from lateral gene transfer events. For ex- 
ample, P. pacificus is the first nematode, which is not a 
plant parasitic species yet contains cellulase genes [25]. P. 
pacificus cellulase genes are most similar to genes found 
in the social amoeba of the Dictyostelid group [45]. These 
cellulases are fully functional and have been maintained 
over long evolutionary time periods [46]. We have identi- 
fied seven cellulase genes in the current transcriptome an- 
notation, six of which are represented on our P. pacificus 
microarrays. In a dauer versus mixed-stage comparison, 
four of these cellulases are significandy differentially 
expressed, all of them being down regulated in dauers (or 
equivalently, expressed predominantly in non-dauer mix- 
stages, Additional file 4: Table S3). This is in agreement 
with their hypothesized function of enabling feeding on 
novel food sources [46]. For example, one gene (Con- 
tig66-snapTAU.4) is already significandy up-regulated as 
early as 12 hours post induction (Additional file 4: 
Table S3, column logFC(D0/D12)). 

Another interesting gene family, which was most likely 
acquired from beetles, is the diapausins. These genes en- 
code small peptides, which are thought to protect the 
dormant beetle against microbial infections in diapause. 
Diapausins provide antifungal activity by acting as Ca2+ 
channel blockers [47]. Out of five members of this gene 
family in P. pacificus genome, three are represented on 
our microarrays. Two of them are significantly up- 
regulated in dauers compared to dauer exit (Additional 
file 4: Table S4, column logFC(D0/D12)) and might be 
important for antifungal immunity in the dauer stage on 
the beetle, while the third gene does not respond to the 
dauer exit program. This gene is in fact downregulated 
in a dauer versus mix-stage comparison (column logFC 
(DO/MixStage), Additional fde 4: Table S4) suggesting 
potential anti-fungal role in other developmental stages. 
Thus, while previous studies in our lab have demon- 
strated that both gene families (cellulases and diapau- 
sins) originate from lateral gene transfer events, we are 
now able to show members of both gene families are de- 
velopmentally regulated in agreement with their pro- 
posed functions: cellulase genes are implicated in 
feeding and are hence downregulated in dauers, and at 
least two diapausin gene appears to have a role in pro- 
tecting hibernating stages. These hypotheses generated 
from our microarray data can further be investigated at 
a functional level. 

One third of all P. pacificus genes are pioneer genes, 
which do not show any obvious sequence similarity to 



other organisms [26]. On average, pioneer genes are sig- 
nificantly higher expressed in the dauer stage in com- 
parison to the remaining genes (p-value < 1.0E-06; 
Additional file 1: Figure S3). We speculate that pioneer 
genes could be especially important to the evolution of 
the dauer stage in P. pacificus, which is essential for its 
necromenic life style. 

Discussion 

The dauer is an ecologically regulated developmental 
stage observed in many free-living nematode species and 
hence provides an interesting model to investigate how 
ecological adaptations are integrated into developmental 
pathways during evolution [1]. Extensive studies in C. 
elegans have uncovered the genetic regulators of dauer 
formation and comparative studies from other free living 
nematodes have begun to provide insights into evolution 
of dauer regulatory genes and pathways [33,34,48]. P. 
pacificus presents an ideal model for comparing dauer 
larvae of free living nematodes. In this study we provide 
a comprehensive comparison of gene expression in the 
dauer stage and dauer exit (12 hour post induction) of 
the two nematode model systems P. pacificus and C. ele- 
gans using the Agilent microarray platform. While we 
are aware of the available C. elegans expression data on 
dauer development (e.g. [5,6,8-14], we felt that it was ne- 
cessary to generate C. elegans data de novo for two rea- 
sons. First, we wanted to use the same platform as that 
used for the analysis of P. pacificus in order to enhance 
the power of a direct comparison. Second, we wanted to 
benefit from technical advances in custom and long nu- 
cleotide microarrays. Based on cross-species comparison 
of transcriptomes of the ecologically important dauer 
stage, we draw four major conclusions from our studies. 

First, we provide a list of similarities and differences 
between P. pacificus and C. elegans, which shows an un- 
expected level of divergence at transcriptome level, even 
though the dauer stage and recovery process appear to 
be developmentally conserved. While this comparison 
allows several evolutionary conclusions (see below), the 
P. pacificus data set on its own can be used as a starting 
point for a functional analysis of the dauer stage and 
dauer exit. This data set represents an invaluable re- 
source given the importance of the dauer stage for the 
ecology of this nematode. The association of P. pacificus 
with scarab beetles is restricted to the dauer stage as 
long as the beetle is alive [27,28]. Only after the beetle's 
death, the nematode exits from the dauer stage to feed 
on the microbes, which develop on the carcass of the in- 
sect. Thus, the P. pacificus dauer stage has a well- 
defined ecological niche and the expression profiles 
described in this study will serve as an entry point to fu- 
ture functional studies. For example, our transcriptomic 
data identifies many P. pacificus specific genes as 
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upregulated in the dauer stage, which implies a potential 
function in adaptations enabling survival on beetles. 

Second, we show that metabolic differences exist be- 
tween both species, with different patterns of regulation 
of genes involved in the central carbon metabolism 
pathways (glycolysis, tricarboxylic acid cycle and oxida- 
tive phosphorylation) of P. pacificus upon dauer exit. 
This difference does not result in obvious phenotypic 
changes during dauer exit, but may be linked to the dif- 
ferences in life history traits. This is especially relevant 
given that P. pacificus is adapted for longevity [43]. 
Under experimental conditions, P. pacificus survives for 
up to one year in the dauer stage, whereas C. elegans N2 
dauer larvae die after approximately 22 weeks. Further 
studies will reveal how much of the observed metabolic 
differences are explained by different life-history traits 
versus differences in rate of development. 

Third, this study provides the first nematode evo-devo 
comparison looking at the downstream consequences of 
homologous developmental processes between species 
belonging to different nematode genera. While detailed 
studies between nematodes as distinct as C. elegans and 
P. pacificus have investigated the regulation of vulva 
and gonad development, sex determination and dauer 
formation [30,34,49,50], most of these studies are con- 
cerned with the regulatory mechanisms rather than the 
"executional programs" of the corresponding develop- 
mental processes. The divergence in the expression pro- 
files of C. elegans and P. pacificus adds an important 
new finding to the growing literature of evo-devo. Previ- 
ous studies have indicated the limited conservation in 
the genetic and molecular control of developmental 
processes in P. pacificus and C. elegans. For example, 
vulva induction relies on different signaling pathways, 
requires a novel regulatory linkage and the acquisition 
of novel protein domains in P. pacificus Wnt signaling 
[32]. This type of result has been discussed as an ex- 
ample for the theory of developmental systems drift, 
which proposes that conserved developmental and mor- 
phological structures can be regulated by largely diverse 
regulatory mechanisms [51]. Considering that gene 
regulatory networks are hierarchically structured, with 
possibly different rates of evolution at the top level 
regulatory genes and the most downstream level of ef- 
fector genes [52], it can be argued that in principle, de- 
velopmental systems may diverge due to differences/ 
drift at any of these levels. Unfortunately, unlike the 
top-level regulatory network, the downstream effector 
programs of vulva development have largely escaped 
identification by developmental genetic approaches and 
have not been easily accessible to transcriptome studies 
as they are single-cell or small group of cell responses. 
Our study circumvents this limitation because dauer 



formation is a "whole body response" of the organism 
to harsh environmental conditions. 

The dauer context is also interesting because in P. 
pacificus, the key transcription factors of the dauer regu- 
latory network, DAF-16 and DAF-12, are conserved 
[33,34]. However, this in itself does not indicate the ex- 
tent to which the downstream targets of the regulatory 
network are subject to evolutionary change. Herein, we 
could demonstrate for the first time that the core down- 
stream execution program of a developmental stage can 
differ tremendously between P. pacificus and C. elegans, 
in spite of conservation of upstream regulators like 
DAF-16 and DAF-12. Thus, these observations make a 
case for extending the concept of developmental systems 
drift to the downstream molecular execution of specific 
developmental stages. 

The fourth conclusion is also related to evolutionary 
theory. While future studies will have to reveal how 
much of the observed differences between P. pacificus 
and C. elegans is really of functional importance, it has 
often been assumed that such differences might simply 
be neutral [23,53]. Gene expression, which is neither 
strongly deleterious nor advantageous, previously termed 
"gratuitous expression" [53], will not be under selection 
and will be free to evolve by drift. Consequently, such 
expression is probably not functional. However, such 
arguments may not apply to the dauer stage since nema- 
todes live off their internal limited energy resources and 
any random or neutral transcriptional activity would di- 
minish these limited resources. 

Comparative studies in developmental genetics have 
driven the studies on evolution of developmental 
mechanisms, and with whole genome sequencing of 
many animal species has now highlighted new facets of 
evolutionary dynamics through comparative genomic 
studies [21]. Since transcriptional regulation is a key 
building block in the genotype to phenotype translation, 
comparative transcriptomic studies add another dimen- 
sion to the analysis of evolutionary processes [23]. Our 
work contributes to the growing set of results from com- 
parative transcriptomics in diverse developmental sys- 
tems (e.g. [54-58]). These studies together span a range 
of conclusions, from high transcriptomic conservation at 
one end, to relatively low conservation in others, sug- 
gesting inherent constraints as well as flexibility in the 
evolution of gene regulatory networks [23]. Future stud- 
ies comparing transcriptomes of homologous biological 
processes in related species, will be important for under- 
standing the role of transcriptome evolution in generat- 
ing animal diversity. Ultimately, this will also reveal the 
extent to which the conserved or divergent expression 
changes are subject to adaptive and non-adaptive forces 
during evolution. 
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Methods 

Worm strains and culture 

We used wild-type strains of two distinct nematode spe- 
cies in all of our experiments. For Caenorhabditis ele- 
gans, we used the N2 (Bristol) strain. For Pristionchus 
pacificus, we used the RS2333 strain (formerly known as 
PS312). For mixed stage cultures, 10 to 15 early adults 
were spotted on 10 cm NGM plates and allowed to grow 
at 20°C for 5 days, and washed off with M9 for RNA 
extractions. The dauers for both species were obtained 
from liquid cultures grown at 25°C. For this, 15 to 20 
mixed stage plates (see above) were washed off and sus- 
pended in a final volume of 500 ml of S-medium in a 
3000 ml Erlenmayer flask. 200ul Nystatin (50 mg/ml 
stock in DMSO) and 500ul Kanamycin (20 mg/ml stock 
solution) was added to prevent fungal and bacterial con- 
tamination. On day 1, 5 and 8 of culture, 12.5 ml of 
OP50 (20%w/v stock in S-medium) was added as food 
source. The cultures were grown for 12 to 14 days in a 
shaker incubator @ 220 rpm, 25°C. Dauers were purified 
from liquid cultures between days 12 to 14. The culture 
was centrifuged to obtain a worm pellet, which was then 
incubated with 1% SDS for 30 minutes to kill any non- 
dauers. Live worms were separated by sucrose flotation, 
and a subsequent precipitation in 15% Ficoll400 resulted 
in a pure dauer pellet (confirmed under a stereo micro- 
scope). The washed dauers were incubated overnight in 
0.1 M NaCl at 20°C, to let them recover from any pos- 
sible stress induced during the harsh dauer purification 
process. The dauers were then collected by centrifuga- 
tion, and used for RNA extraction and for starting 
dauer-exit cultures. 

Dauer exit time course 

For the dauer-exit 12-hour time-point samples, around 
200 purified dauers were spotted on each 10 cm NGM 
plate pre-spotted with 2 ml OP50, and grown at 20°C for 
12 hours. These conditions ensured abundant food supply 
and non-crowded conditions, favorable for inducing dauer 
exit [3]. Worms from 20 to 30 such plates were used for 
RNA extraction per sample. The 12-hour time-point was 
chosen for both the species to give them sufficient time to 
recover from dauer stage and manifest measurable tran- 
scriptional changes. The recovering worms from both spe- 
cies were monitored to ensure comparable developmental 
rates based on the following criteria: (1) Most of the popu- 
lation (~90%) of the recovering worms from C. elegans as 
well as P. pacificus resume pharyngeal pumping within 3 
hours after inducing dauer exit. (2) At the 12 hour time- 
point (the stage of dauer-exit samples), no discernible mor- 
phological differences are found between the recovered 
dauers from both species. (3) The recovering animals from 
both the species enter a lethargus between 13 to 14 hours 
post-recovery [3]. (4) The worms that are allowed to 



develop further at 20° C reach the next moult (L4 for 
C. elegans, J4 for P. pacificus) between 22 to 23 hours after 
induction of dauer exit. (5) Finally, recovered worms from 
both species started laying eggs between 42 to 45 hours 
post-induction. 

RNA extraction 

Total RNA was isolated using TRIzol (Invitrogen) 
according to the manufacturer's instructions. Four bio- 
logical replicates for each stage (mixed, dauers, dauer- 
exit at 12 hrs) were prepared from parallel cultures. The 
extracted RNA was purified by using phenolxhloroform: 
isoamyl alcohol precipitation, to remove any trace of 
TRIzol contamination which might interfere with subse- 
quent reactions. The RNA pellet was suspended in 
RNAse free water, quantity and quality was assayed on a 
Nanodrop spectrophotometer, and stored at -80°C until 
the microarray reactions. 

Microarray design and experiments 

We used the Agilent C. elegans oligonucleotide microar- 
rays (~ 43,000 probes for ~ 20,000 open reading frames, 
GEO accession: GPL10094) for all C. elegans experi- 
ments. For P. pacificus, we designed custom Agilent 
microarrays based on the most recent transcriptome an- 
notation (~ 23,000 gene predictions, predicted with ex- 
ternal evidence from 454 EST alignments). We could 
accommodate ~ 93,000 probe sequences on our custom 
P. pacificus microarrays (NCBI GEO platform accession 
GPL14372). This probe set was designed using the Oli- 
goWiz 2.0 software [59] and submitted for custom fabri- 
cation to Agilent Technologies via their eArray web tool. 
Out of 93,000 probes, 87,070 probes map uniquely to 
the latest genome assembly [26]. We restrict all subse- 
quent analyses to a probe set of 69,916 high-confidence 
probes. 

Microarray hybridizations were carried out in a two- 
color format, with 4 biological replicates per experiment 
including a pair of dye-flips. Two sets of hybridizations 
were carried out: (1) Dauers versus Mix-Stage, and (2) 
Dauer-Exit at 12 hr versus Mix-stage, thus making a total 
of 8 microarray hybridizations per species. Equal amounts 
of total RNA (500 ng to 1000 ng) for 4 biological replicates 
from each stage was used to prepare Cy5- or Cy3- labelled 
cRNA using Quick Amp Labelling Kit (Agilent Technolo- 
gies, Inc, Santa Clara, CA, USA), as per manufacturer's 
instructions. Based on the starting amounts of total RNA, 
appropriate dilutions of positive controls (Spike Mix-A and 
Spike-Mix B from Agilent) had been added to the 
reaction-mix before the RT reaction, as per manufacturer 
instructions. The labelled cRNAs were hybridized either to 
Agilent C. elegans oligonucleotide microarrays or our cus- 
tom P. pacificus microarrays from Agilent. The arrays were 
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scanned using GenePix 4000B Microarray Scanner, and 
raw data extracted using GenePix Pro sofware (version 6). 

Analysis of microarray expression data 

The data was analyzed using the Bioconductor [60] 
package limma [61]. Briefly, the array quality was 
assessed by checking for uniform background and fore- 
ground intensities over the entire array. The signal was 
background corrected using the normexp method [62]. 
The arrays were lowess normalized individually ("nor- 
malizeWithinArrays" option), with differential weights 
assigned to probes and to positive control spike-ins, 
which are expected to show no fold change [63]. This 
differential weighing of probes is particularly necessary 
to account for changes in the relative proportion of 
mRNA versus total RNA. Without this differential 
weighing scheme, the fold change calculations can be er- 
roneous [64]. The weight parameters were optimized 
based on MA-plots such that spike-in controls show 
their expected fold change values. lmFit function was 
used to fit a linear model to probe intensities across 
arrays, and differential expression was calculated by em- 
pirical Bayes method using the eBayes function [65]. 
Control of FDR was employed as correction for multiple 
testing. All the data from this publication have been 
deposited in a MIAME compliant format [66] at NCBIs 
Gene Expression Omnibus database (http://www.ncbi. 
nlm.nih.gov/geo/) and are accessible through GEO Series 
accession numbers GSE30977 and GSE31861. 

Mapping of 1:1 orthology relations 

We used a pairwise best BLASTP strategy to compute 
1:1 orthologs. Briefly, we ran all protein sequences from 
C. elegans as query versus the database of P. pacificus 
gene predictions and vice versa. Only hits with a BLAST 
score > = 50 bits were retained. We define mutually best 
hits as 1:1 orthologs. We identified 7,176 ortholog pairs 
with this methodology. 

Gene ontology enrichment analysis 

GO ontologies for C. elegans were downloaded from 
wormmart WS200. Gene ontologies were assigned to P. 
pacificus genes based on the previously defined 1:1 
orthologs. The topGO tool [37] (version 2.4.0) was used 
for computing significantly enriched GO terms. We used 
the "GOFisher" test statistic and a p-value cutoff of 0.05. 
The background set was limited to the 6,126 orthologs 
that are represented on both C. elegans and P. pacificus 
microarrays. We did not apply any multiple testing cor- 
rection to the reported p-values. We followed this strat- 
egy to uncover all "trends" in the data to attain a 
comprehensive picture of the underlying biology. 



Expression cluster enrichment analysis 

WormBase [36] contains information on co-expressed 
gene groups in C. elegans. The list of microarray experi- 
ments where a given C. elegans gene is known to be dif- 
ferentially expressed can be extracted from the section 
"Expression Cluster" from the WormBase gene summary 
page for each gene. We retrieved all available expression 
clusters for C. elegans genes from the WormBase web 
site. We inferred expression clusters for P. pacificus 
based on the set of 1:1 orthologs. Only the clusters asso- 
ciated with the 6,126 1:1 orthologs were used as the 
background set. P-values for expression cluster enrich- 
ment in dauer-enriched or dauer exit-enriched gene sets 
were computed with a 2x2 Fisher exact test, and a mul- 
tiple testing correction to control FDR was applied. 

Pfam domain annotation 

We annotated the proteome of C. elegans and P. pacifi- 
cus with PFAM domain matches (PFAM release V23/4 
[67])). Protein HMM searches were performed with 
HMMer 3.0 [68] using a p-value cutoff of 0.001. 

KEGG pathway analysis 

We retrieved the latest gene to pathway mapping for C. 
elegans from the KEGG SOAP server (see www.genome. 
jp/kegg/soap for details). Again, pathway annotations 
were transferred to P. pacificus with the help of the 1:1 
orthologs. Gene expression information and pathway in- 
formation were overlaid to generate Figure 3. 

Data access 

Microarray data have been deposited in NCBIs Gene Ex- 
pression Omnibus database (http://www.ncbi.nlm.nih. 
gov/geo/) and are accessible through GEO Series acces- 
sion numbers GSE30977 and GSE31861. The gene set 
annotation for this manuscript may be obtained from 
http://www.pristionchus.org/download/ . The list of 1:1 
orthologs and PFAM domain annotations may be also 
obtained from there. 

Additional files 



Additional file 1: Figure SI. The impact of different normalization 
strategies on expression fold change estimation, a) The general 
concept of using spike-ins to control for changes in total mRNA 
abundance relative to total RNA. See van de Peppel ef al., 2003 for a 
detailed explanation, b) MA-plot of a Loess based normalization 
assuming a net mRNA expression fold change of 0 between the two 
samples. The M axis depicts the log2 fold change between the red and 
green microarray channels. The A axis is the average log2 intensity across 
both microarray channels. The position of the spike-in signals {colored 
dots) clearly indicate a deviation from this assumption. The 
corresponding legend in the upper left shows the expected fold change 
of the spike-ins relative to total RNA levels. Consequently, the location of 
the spike-ins corresponds to a down-shift of mRNA abundance levels 
within the total RNA population, c) Differentially expressed gene counts 
for the dauer vs. mixed stage comparison conditional on the 
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normalization strategy. Figure S2. Explanation of a common reference 
design for microarrays (see Eisen and Brown 1999). For each species, RNA 
from biological replicates of mix-stage worms (M1 to M4) were 
combined together to generate one common reference pool. Labelled 
aRNA produced from independent biological replicates were then co- 
hybridized with labelled aRNA from the common reference pool, 
including two dye-swaps. Samples D1 to D4 represent Dauer samples 
and DEI to DE4 represent dauer exit samples at 12 hour timepoint. The 
blue arrows indicate the direction of labelling in each co-hybridization 
(arrow head = Cy3 labelled, arrow tail = Cy5 labelled). This design allows 
comparison via the common reference pool, but at the same time 
remains flexible for adding more time-points to the study, if needed. 
Figure S3. Cumulative plot of log2 expression fold changes for genes 
with and without sequence conservation (pioneer genes). The two 
distributions are significantly different (two-sample Kolmogorov-Smirnov 
test; p-value<10E- 16). 

Additional file 2: Table SI. Differential expressed genes in the 
dauer versus dauer-exit comparison in P. pacificus. 

Additional file 3: Table S2, Differential expressed genes in the 
dauer versus dauer-exit comparison in C. elegans. 

Additional file 4: Table S3. Differential expression of P. pacificus 
cellulase genes that were acquired by horizontal gene transfer. 

Table S4. Differential expression of P. pacificus diapausin genes that were 
acquired by horizontal gene transfer. 
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